classdef EddyArrayMatchSSTAClass < handle
    % 需要初始化温度的读取地址
    properties
        eddyArray
        % type 没用上
        type
        extent = [0 25;108 122];
        year = [1993 2020];
        sstaFilePath = 'F:\按天气候态\'
    end
    
    methods
        function obj = EddyArrayMatchSSTAClass(eddyArray, varargin)
            obj.eddyArray = eddyArray;
            idx_sstaFilePath = find(strcmp(varargin,'sstaFilePath'));
            if ~isempty(idx_sstaFilePath)
                obj.sstaFilePath= varargin{idx_sstaFilePath + 1};
            end            
        end

        % 计算区域EKE
        function [lon,lat,EKE] = area_EKE_measure(obj)
            beginTensor = datetime(obj.year(1),1,1);
            endTensor = datetime(obj.year(2),12,31);
            t = beginTensor;
            EKE_sum = 0;
            EKE_count = 0;
            while(t<=endTensor)
                eddy.time = t;
                [lon, lat,u, v, ~] = EddyArrayMatchSSTAClass.u_v_select_by_eddy(eddy);
                lonidx = lon > EAMSC.extent(2,1) & lon < EAMSC.extent(2,2);
                latidx = lat > EAMSC.extent(1,1) & lat < EAMSC.extent(1,2);
                lon = lon(lonidx);
                lat = lat(latidx);
                u = u(lonidx,latidx);
                v = v(lonidx,latidx);

                EKEItem = (u.^2 + v.^2)/2;
                EKE_sum = EKE_sum + EKEItem;
                EKE_count = EKE_count +1
                t = t + days(1);
            end
            EKE = EKE_sum./EKE_count;      
        end
        
        % 计算 涡致EKE/总EKE 表征涡旋涡动能贡献率
        function [lon,lat, LEP] = area_LEP_measure(obj)
            beginTensor = datetime(obj.year(1),1,1);
            endTensor = datetime(obj.year(2),12,31);
            t = beginTensor;
            EKE_sum = 0;
            EKE_count = 0;
            while(t<=endTensor)
                eddys = obj.eddyArray(find([obj.eddyArray.time] == t));
                [lon, lat,u, v, ~] = EddyArrayMatchSSTAClass.u_v_select_by_eddy(eddys(1));
                lonidx = lon > obj.extent(2,1) & lon < obj.extent(2,2);
                latidx = lat > obj.extent(1,1) & lat < obj.extent(1,2);
                lon = lon(lonidx);
                lat = lat(latidx);
                u = u(lonidx,latidx);
                v = v(lonidx,latidx);
                [m,n] = size(u);
                mask = zeros(m,n);

                for i = 1:length(eddys)
                    [x,y] = meshgrid(lon, lat);
                    x = x';
                    y = y';
                    longitude = eddys(i).contour_longitude;
                    latitude = eddys(i).contour_latitude;
                    in = inpolygon(x,y,longitude,latitude);
                    mask(in==1) = 1;
                    % longitude = eddys(i).longitude;
                    % latitude = eddys(i).latitude;
                    % for p = 1:m
                    %     for q = 1:n
                    %         in = inpolygon(lon(p),lat(q),longitude,latitude);
                    %         if in == 1
                    %             mask(p,q) = 1;
                    %         end
                    %     end
                    % end
                end
                EKEItem = 0.5*(u.^2 + v.^2).*mask;
                EKE_sum = EKE_sum + EKEItem;
                EKE_count = EKE_count +1
                t = t + days(1);
            end
            LEP = EKE_sum./EKE_count;      
        end

        
        function [SSTA, DhA, R, depth] = mean_eddyMatchedSSTA_Radis(obj)
            eddyMatchedSSTA_Sum = 0;
            eddyMatchedSSTA_Count = 0;
            eddyMatchedSSTA_Radis_Square_Sum = 0;
            eddyMatchedDhA_Sum = 0;
            eddyMatchedDhA_Count = 0;
            depth = nan;
            for eddyArrayi = 1:length(obj.eddyArray)

                disp(['已经计算' num2str(eddyArrayi) '个涡旋的平均']);
                eddy = obj.eddyArray(eddyArrayi);
                try
                    [eddyMatchedSSTA, eddyMatchedDhA, Radis, depth] = EddyArrayMatchSSTAClass.eddyMatchedSSTA_measure(obj.sstaFilePath, eddy);
                catch
                    % warning(['第' num2str(eddyArrayi) '个涡旋匹配不到']);
                    continue;
                end

                %%% 添加限制条件
                if EddyThermalTransortClass.thermalMeasure(eddyMatchedSSTA(21:41,21:41,:), Radis, depth) * eddy.type > 0
                    continue;
                end
                %%%

                eddyMatchedSSTA_Radis_Square_Sum = eddyMatchedSSTA_Radis_Square_Sum + Radis^2;
                [eddyMatchedSSTA_Sum, eddyMatchedSSTA_Count] = sum_measure(eddyMatchedSSTA_Sum, eddyMatchedSSTA_Count, eddyMatchedSSTA);
                [eddyMatchedDhA_Sum, eddyMatchedDhA_Count] = sum_measure(eddyMatchedDhA_Sum, eddyMatchedDhA_Count, eddyMatchedDhA);


            end

            eddyMatchedSSTA_mean = eddyMatchedSSTA_Sum ./ eddyMatchedSSTA_Count;
            eddyMatchedDhA_mean = eddyMatchedDhA_Sum ./ eddyMatchedDhA_Count;
            eddyMatchedSSTA_Radis_mean = sqrt(eddyMatchedSSTA_Radis_Square_Sum / length(obj.eddyArray));
            DhA = eddyMatchedDhA_mean;
            SSTA = eddyMatchedSSTA_mean;
            R = eddyMatchedSSTA_Radis_mean;

            function [Sum, Count] = sum_measure(Sum, Count, data)
                Count_Matrix = ones(size(data));
                nanIndex = find(isnan(data));
                data(nanIndex) = 0;
                Count_Matrix(nanIndex) = 0;
                Count = Count + Count_Matrix;
                Sum = Sum + data;
            end
        end
        % 计算涡旋EKE
        function [EKE] = mean_eddyMatchedUV_EKE(obj)
            EKE_Sum = 0;
            EKE_Count = 0;
            for i = 1:length(obj.eddyArray)
                eddy = obj.eddyArray(i);
                try
                    [eddyMatchedU, eddyMatchedV] = EddyArrayMatchSSTAClass.eddyMatchedUV_measure(eddy);
                    EddyMatchedEKE = (eddyMatchedU.^2 +eddyMatchedV.^2)/2*10000;%m转cm
                    meanEKE = mean(mean(EddyMatchedEKE));
                    disp(['已经计算' num2str(i) '个涡旋的平均']);
                catch
                    continue;
                end
                EKE_Sum = EKE_Sum + meanEKE;
                EKE_Count = EKE_Count + 1;
            end
            EKE = EKE_Sum/EKE_Count;
        end
    end
    methods(Static)
        % 涡旋与三维温盐场匹配工作
        function [eddyMatchedSSTA, eddyMatchedDhA, Radis, depth] = eddyMatchedSSTA_measure(filepath, eddy)
            [lon, lat ,depth ,ssta, dha] = EddyArrayMatchSSTAClass.ssta_select_by_eddy(filepath, eddy);
            Radis = eddy.radius;
            % 地球半径
            R = 6371000;
            
            lonDeltaDistance = pi/180*R;
            % latDeltaDistance = 111000 * cosd(eddy.latitude);    
            delta = Radis/lonDeltaDistance;

            interpx = -3:0.1:3;
            interpy = -3:0.1:3;
            
            interplon = double(interpx * delta + eddy.longitude);
            interplat = double(interpy * delta + eddy.latitude);
            [x,y,z] = meshgrid(lat, lon, depth);
            [xx, yy, zz] = meshgrid(interplat, interplon, depth);
            eddyMatchedSSTA = interp3(x,y,z,ssta,xx,yy,zz);
            eddyMatchedDhA = interp3(x,y,z,dha,xx,yy,zz);
        end

           
        function [lon, lat ,depth ,ssta, dha] = ssta_select_by_eddy(filepath, eddy)
            eddy_date = eddy.time;
            MatchSSTFullFileName = selectMatchSSTFullFileName(filepath, eddy_date);
            load(MatchSSTFullFileName,'lon','lat','depth','ssta','dha');
                        % 根据时间读取ssta

            function MatchSSTFullFileName = selectMatchSSTFullFileName(filepath, eddy_date)
                year = eddy_date.Year;
                % MatchSSTFullFileName = ['G:\GOPR按月气候态ssta\' num2str(year) '\' datestr(eddy_date, 'yyyymmdd') '.mat'];
                MatchSSTFullFileName = [filepath num2str(year) '\' datestr(eddy_date, 'yyyymmdd') '.mat'];
            end

            % cemes周数据
            % function MatchSSTFullFileName = selectMatchSSTFullFileName(eddy_date)
            %     year = eddy_date.Year;
            %     armorfilePath = 'H:\code\大论文数据处理\数据\armor3D_new\armor3D季节平均位势温度异常计算\beforefilter\';
            %     sst_file_list = [dir([armorfilePath num2str(year - 1) '\' '*.mat']);
            %                     dir([armorfilePath num2str(year) '\' '*.mat']);
            %                     dir([armorfilePath num2str(year + 1) '\' '*.mat'])];
            % 
            %     for i = 1:length(sst_file_list)
            %         fileName = sst_file_list(i).name;
            %         date_str = fileName(1:8);
            %         sst_date = datetime(date_str,'InputFormat','yyyyMMdd');
            %         if abs(eddy_date - sst_date) <= 3
            %             MatchSSTFullFileName = [sst_file_list(i).folder '\' sst_file_list(i).name];
            %             break;
            %         end
            %     end
            % end
        end

        function [eddyMatchedU, eddyMatchedV] = eddyMatchedUV_measure(eddy)
            [lon, lat ,u, v, ~] = EddyArrayMatchSSTAClass.u_v_select_by_eddy(eddy);
            Radis = eddy.radius;
            % 地球半径
            R = 6371000;
            
            lonDeltaDistance = pi/180*R;
            % latDeltaDistance = 111000 * cosd(eddy.latitude);    
            delta = Radis/lonDeltaDistance;

            interpx = -1:0.1:1;
            interpy = -1:0.1:1;
            
            interplon = double(interpx * delta + eddy.longitude);
            interplat = double(interpy * delta + eddy.latitude);
            [x,y] = meshgrid(lat, lon);
            [xx, yy] = meshgrid(interplat, interplon);
            eddyMatchedU = interp2(x,y,u,xx,yy);
            eddyMatchedV = interp2(x,y,v,xx,yy);
        end

        function [lon, lat,u, v, sla] = u_v_select_by_eddy(eddy)
            eddy_date = eddy.time;
            MatchUVFullFileName = selectMatchUVFullFileName(eddy_date);
            load(MatchUVFullFileName,'longitude','latitude','ugosa','vgosa','sla');
            lon = longitude;
            lat = latitude;
            u = ugosa;
            v = vgosa;
            
            function MatchUVFullFileName = selectMatchUVFullFileName(eddy_date)
                year = eddy_date.Year;
                MatchUVFullFileName = ['H:\code\南海\SLA_ok\' num2str(year) '\' datestr(eddy_date, 'yyyymmdd') '.mat'];
            end
        end
    end
end

